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SUBROUTINE  PROBDIF 


INTRODUCTION 


Given  n  independent  points  ,  X2 ,  Xj^  from  the  density  ')  and  m 

Independent  points  y^,  y2,  y^i  from  the  density  we  wish  to  deter- 

mine  probability  that 


f^(Z)  -  g(Z) 
fy(Z)  «  g(Z-ij)  . 

That  is,  we  wish  to  determine  the  probability  chat  one  density  is  shifted  by 
U  from  the  ocher  density.  Since  we  do  not  want  to  assume  anything  about  the 
density  g{*),  we  will  employ  a  nonparamecric  procedure. 

GENERAL  DESCRIPTION  OF  TECHNIQUE 

If  the  x's  and  y's  come  from  the  same  density,  there  are  (m+n)!/n!m! 
possible  orderings  of  the  points;  and  each  ordering  is  equally  likely.  If 
there  are  n'  x's  smaller  than  any  y  and  m'  y's  larger  than  any  x,  the  prob¬ 
ability  Chat  the  x's  and  y's  come  from  a  different  density  is  equal  to  the 
probability  of  obtaining  any  sequence  that  starts  with  n '  or  more  x’s  and 
ends  with  m"  or  more  y's. 

As  an  example  of  this  procedure,  let  us  assume  that  there  are  3  x  sam¬ 
ples  1.0,  2.0,  and  3.0  and  3  y  samples  3.5,  ■i.S,  and  5.5.  Given  these 
examples,  there  are  20  possible  orderings  for  the  samples  and  these  order¬ 
ings  are  listed  in  Table  1.  The  probability  chat  the  x  points  are  u  less 
than  Che  y  points  is  given  in  Table  2.  For  u-0  the  sequence  has  3  x's 
followed  by  3  y's.  Only  sequence  IH  satisfies  this  condition  and  its 
probability  of  occurrence  is  0.05.  Thus,  Che  probability  that  the  samples 
do  not  come  from  the  same  density  is  one  minus  this  probability  or  0.95. 

For  1^1,  the  sequence  starts  with  2  x's  and  ends  with  2  y's.  Sequences 
’>1  and  ?/2  satisfy  this  condition,  their  probability  of  occurrence  is  9.10, 
and  the  probability  that  Che  samples  do  not  come  from  the  same  density  is 
0.90.  For  'J»2,  the  sequence  starts  with  an  x  and  ends  with  a  y.  Sequences 
'^1 ,  d2,  ’^3,  •■^5,  9b,  and  ‘g  satisfy  this  condition,  their  probability  is 
0.3,  and  thus  Che  probability  that  the  samples  do  not  come  from  the  same 
density  is  0.7.  When  the  sequence  starts  with  y's  and  ends  in  x's,  the 
probability  of  occurrence  of  the  sequences  equals  the  probability  chat  the 
samples  do  not  come  from  Che  same  density;  i.e.,  one  does  not  have  to  sub¬ 
tract  the  sequence  probabilities  from  1. 

When  Che  number  of  samples  is  a  moderate  number,  it  becomes  irapraccicai 
CO  enumerate  all  possible  sequences.  For  instance,  if  n*m“15,  Che  number  of 
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Table  1  -  Possible  Sequences 


Table  2 


Sequence 

Sequence  I 

1 

1 

xxxyyy 

2 

xxyxyy 

3 

xxyyxy 

4 

xxyyyx 

5 

xyxxyy 

6 

xyxyxy 

7 

xyxyyx 

3 

xyyxxy  1 

9 

xyyxyx  j 

10 

xyyyxx 

11 

yxxxyy  ! 

12 

yxxyxy 

13 

yxxyyx 

14 

yxyxxy 

15 

yxyxyx 

16 

yxyyxx 

17 

yyxxxy 

18 

yyxxyx 

19 

yyxyxx 

yyyxxx 

-  Probabillcy  chat  Che  density  of  Che  x  points  shifted  by 
U  is  less  than  the  density  of  the  y  points. 


u 

!  1 

j  Probability  1 

i 

1  0 

i  1 

1  -55  , 

i  i 

.90  j 

1  2 

1  .70  ' 

j  3 

j  .30  j 

1 

.10  ; 

1  5 

t 

j  .05  i 

1  1 

possible  sequences  is  greater  than  10®.  Consequently,  we  will  discuss  a 
different  technique  for  calculating  Che  probability  of  the  various  sequences 
without  enumerating  them.  This  technique  will  also  handle  the  cases  where 
Che  smallest  and  largest  samples  come  from  the  same  sec  of  points. 

SPECIFIC  TECHNIQUE 

Let  us  first  define  as  the  number  of  x's  which  are  smaller  Chan  all 
of  the  y's  and  NSy  as  the  number  of  y's  which  are  smaller  chan  all  of  the 
x's.  Similarly,  define  as  Che  number  of  x's  which  are  larger  than  all 
of  the  y's  and  NLy  as  the  number  of  y's  which  are  larger  than  all  of  the 
x's.  Note  that  if  NSjt  >  0,  then  NSy  »  0;  and  if  NSy  >  0,  then  NS^  =  0. 
Similarly,  if  NL^  >  0,  then  NLy  “  0;  and  if  NLy  >  0,  then  NL^  =  0.  Based  on 
these  four  values  (NS^,  NSy,  NL^  and  NLy),  four  different  cases  are  possible 

1.  NSx  >  0  and  NLy  >  0, 

2.  NSy  >  0  and  NL^  >  0, 

3.  NSx  >  0  and  NLx  >  0. 

4.  NSy  >  0  and  NLy  >  0. 

Each  of  these  cases  will  now  be  discussed  and  a  formula  will  be  generated 
to  calculate  the  probability  that  the  x  data  set  comes  from  a  density  which 
has  a  smaller  location  parameter  than  the  y  density. 

Case  1 ;  NSy  >  0  and  NLy  >  0 

If  NS^  •  n,  that  la  all  the  x's  are  less  than  all  the  y’s,  the  proba¬ 
bility  of  obtaining  Che  first  x  in  the  sequence  is  n/(m+n),  the  probability 
of  obtaining  the  second  x  in  the  sequence  is  (n-l)/m+n-l),  and  the  proba¬ 
bility  of  Che  third  x  is  (n-2)/m+n-2) ,  and  the  probability  of  the  i-th  x  is 
<n+l-i)/(m+n+l-i) .  Consequently,  the  probability  of  this  sequence  is 

n  .  n-1  .  n-2  , , ,  1  «  n!m! 

m+n  m+n-1  at+n-2  m+l  (m+n)!  (1) 

If  NS^  <  n,  let  1  be  in  Che  range  NS^  ^  i  <  n;  and  let  i  be  the  exact 
number  of  x's  which  start  a  sequence.  The  probability  that  exactly  i  x’s 
start  a  sequence  follows  from  the  above  argument  and  equals 

n  ,  n-1  . , .  n-t-l-i  ,  n!  (m-t-n-i) ! 

m+n  m+n-1  m+n+l-i  (n-i)!  (m+n)!  (2) 

Since  Che  sequence  starts  with  exactly  i  x's,  Che  next  sample  must  be  a  y; 
and  Che  probability  of  this  y  is 

m 

(m+n-1)  (3) 


So  far  we  have  considered  1+1  samples;  1  samples  coming  from  the  x  sec  and 
1  sample  coming  from  the  y  set.  Let  us  now  assume  Chat  the  sequence  ends 
with  exactly  j  y  samples.  This  sequence  of  y  samples  must  be  proceeded  by 
an  X,  and  the  probability  of  this  x  is 


n-1 

nH-n-i-1 

Finally  the  probability  of  ending  with  j  y  samples  is 

m-1  .  m-2  , ,  ^  m-.i  .  (arl) !  (iiri-n-i-j-2)  ’■ 

ittfn-i-2  iiH-n-i-3  m+n-l-j-1  (m-j-l)!  (iii+n-i-2)! 


Thus,  the  probability  P^j  of  all. the  sequences 
and  end  with  exactly  j  y's  is  the  product  of  e 
(5) ,  i.e. , 


s  that  start  with  exactly  i  x's 
expressions  (2),  (3),  (4),  and 


p  «  nt  ( m+n-i ) !  m  n-1  (m-l) !  (iiH-n-i-.j-2) ! 

(n-1)!  (m+n)!  (aH-n-l)  (nrt-n-i-l)  (m-j-l)!  (iiH-n-i-2)!  (6) 

Simplifying  this  expression,  one  obtains 

p  s  n!m!  1 _  (iiH-n-l-j-2) ! 

(m+n)!  (n-1-1)!  (m-j-l)!  (7) 

For  example,  if  m^n--3,  1”2,  and  j“l,  ■  1/20,  which  corresponds  to 

sequence  #3  in  Table  1.  If  m-n*3,  !■!,  and  j“l;  P^j  »  2/20,  which  corre¬ 
sponds  to  sequences  //6  and  #8  in  Table  1.  Equation  (7)  is  valid  for  NS^  £ 
i  <  n  and  NLy  ^  j  <  m.  Thus,  the  probability  of  obtaining  sequences  which 
start  with  at  least  NS^  x's  and  end  with  at  least  NLy  y's  can  be  obtained  by 
sxunming  Eq.  (7)  over  the  proper  ranges  of  i  and  j  and  adding  expression  (1) 
to  this  sum.  Consequently,  by  subtracting  this  probability  from  1,  one 
obtains  the  probability  that  the  x  points  come  from  a  density  which  has  a 
smaller  location  parameter  that  the  y  density  has.  The  desired  probability 
is 


P  -  1  - 


(m+n) ! 


( n-1-1) ! 


(m+n-i-j-2) ! 
( m- j -1 ) ! 


Case  2 ;  NSv  >  0  and  ML*  >  0 


This  case  is  the  dual  of  case  1.  If  we  let  NSy  ^  i  <  m  and  j  <  n, 

the  revised  Eq.  (8)  would  yield  the  probability  that  the  y's  are  less  chan 
the  x's.  The  probability  that  the  x's  are  less  than  the  y's  is  1  minus  this 
probability  or 


( n-Hn) ! 


(m-i-D! 


E 


(m+n-i-j-2)  ! 
(  n- j  +1 )  ! 


Case  3 ;  NS*  >  0  and  Nl~  >  0 


This  case  is  likely  to  occur  if  n  >>  m.  If  the  desired 

probability  is 
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P  =■  0.5 


<  10) 


If  SSx  ^^Lx,  iat  MS  =  Minimum  (NSx.  MLx)  ML  =  Maximum  (cMSx,  MLx)- 

Then,  Che  sequences  of  interest  are  chose  sequences  which  start  with  ML  jr 
more  x's  and  end  with  either  MS  or  less  x’s  or  any  number  of  y's.  The 
probability  chat  a  sequence  starts  with  n  x's  and  ends  with  m  y's  is  given 
by  expression  (1).  Let  i  be  in  the  range  ML  ^  i  <  n  and  let  i  be  the  exact 
number  of  x's  which  start  a  sequence.  The  probability  that  exactly  i  x's 
start  a  sequence  is  given  by  (2).  Since  che  sequence  starts  with  exactly  i 
x's,  the  next  sample  must  be  a  y;  and  Che  probability  of  this  y  is  given 
by  (3).  So  far  we  have  considered  i+1  samples:  i  samples  coming  from  che 
X  set  and  1  sample  coming  from  the  y  set.  Let  us  now  assume  chat  che 
sequence  ends  with  exactly  j  x  samples.  This  sequence  of  x  samples  must 
be  proceeded  by  a  y,  and  the  probability  of  this  y  sample  is 

m+n-1-1  (11) 


The  probability  of  ending  with  j  x  samples  is 

n-i  ,  n-i-1  . . .  n-l-j-*-l  ,  ( n-1 ) !  (m-t-n-i-j-2) ! 

m+n-1-2  m+n-i-3  m+n-i-j-1  (n-i-j)!  (m+n-i-2)!  (12) 


The  possible  values  of  j  are  from  1  to  wjere  M]^  is  che  minimum  of  MS  and 
(n-i).  The  probability  of  starting  a  sequence  with  exactly  i  x's  and  end¬ 
ing  it  with  exactly  j  x's  is  che  product  of  (2),  (3),  (11),  and  (12) 


nl  (m-fn-1) !  m  m-1  (n-i) ! 

(n-i)!  (m+n)!  (m+n-i)  (am\-t-l)  (n-i-j)! 


(m+n-i-j-2) I 
(m+n-i-2) ! 


which  simplifies  to 


n!m(m-l)  ( m+n-i -j -2)! 
(m+n)!  (n-i-j)! 


( 13) 


(14) 


The  sequences  of  interest  which  start  with  exactly  i  x's  can  also  end  with 
j  y  samples  and  this  probability  is  given  by  (7).  Combining  (7)  and  (14) 


and  summing 

over  the 

proper  values  of 

1  and  j ,  one  obtains 

p'- 

n!m! 

n-l  1 

*  2 

i»NL 

f 

1  n!m(m- 

1)  \  '  (m+n-i-j -2)  ! 

(m+n) ! 

1  (oi+n) 

!  /  J  (n-l-j)! 

j=l 

+ 

n!m! 

1 

ti-L 

z 

(m+n-i-i-2)!  | 

1 

( m+n ) ! 

( n-i-i) ! 

(m-j-i)!  1 

r 

Rearranging  terms,  one  obtains 
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i  a-i 

r 

P*  -  1  1  ^ 

(m-t-n) !  j 

(m-2/ !  /  J 

f  i-NL 

j-1 

m-1 

.1..  ..  y* 

(m-m-l-i-2) ! 

(n-i-l)!  ^ 

(m-j-1) ! 

j-1 

. 

n-  i- j  -  2 ) ! 


Finally,  If  >  >0^ 

P  -  1  -  P'  ; 


and  if  KL^  >  NS^ 


P  -  P* 


(16) 


(17) 


18 ) 


Case  4:  MSy  >  0  and  NLy  >  0 

This  case  is  the  dual  of  case  3.  Letting  MS  =•  Minimum  (MSy,  MLy) ,  ML  = 
Maximum  (MSy,  MLy),  and  M2  be  Che  minimum  of  NS  and  (m-i),  the  probability 
P'  is  given  by 


P  - 


(m+n)  i 


m-1  1 

r  ■'^2 

‘  *  r  1 

1  ^ 

(n-2)'! 

i-NL  1 

L  j-1 

n-1 

-] 

1  V 

(m+n-i-j-2) ! 

(m-l-l)' 

(n-j-l)' 

j-1 

J) 


(19) 


If  MSv  >  MLy, 


P  =•  P'  ; 


20) 


and  if  ML,,  >  MS,, 


P  =•  1  -  P'  . 


'll 


SUBROUTINE  PROBDIF 

Depending  on  the  values  of  MS^,  MS,,,  ML^,  and  MLy,  the  desired  proba¬ 
bility  is  given  in  Sq.  (3),  Eq .  (9),  Eqs.  (ib-l3),  or  Eqs.  (19-21).  All 
of  these  equations  involve  douole  summacxuns  over  factorials  anl  the  com¬ 
putation  time  will  vary  as  n^  (or  m^)  if  these  formula  are  implemented 
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directly.  However,  since  the  terms  involve  factorials,  iterative  formulas 
can  be  developed  and  the  computation  time  can  be  made  proportional  to  n-. 
Such  a  technique  was  incorporated  into  subroutine  PROBDIF  which  calculates 
the  desired  probability. 

The  calling  sequence  for  the  subroutine  is 

CALL  PROBDIF  ( nx,x, ny ,y ,xdb, prob, Iset) 


where 


nx  is  number  of  points  from  the  x  distribution 
x(  .)  are  points  from  the  x  distribution 
ny  is  number  of  points  from  the  y  distribution 
y(  .)  are  points  from  the  y  distribution 
xdb  is  the  db  difference  to  be  tested 
prob  is  the  probability  of  a  xdb  difference 
Iset^l  indicates  that  points  have  already  been 
ordered  from  smallest  to  largest. 

Note,  if  successive  calls  are  made  to  PROBDIF  with  different  values  of 
xdb,  iset  should  be  set  to  1  after  the  first  call  to  minimize  computation 
time.  A  Fortran  listing  of  PROBDIF  is  given  in  appendix  A. 

SUMMARY 


'  Given  two  sets  of  data,  this  report  describes  a  technique  for  deter¬ 
mining  the  probability  that  the  first  density  is  a  shifted  version  of  the 
second  density.  We  assumed  nothing  about  the  density  and  used  a  uunpara- 
metric  procedure  based  on  the  smallest  and  largest  samples  of  the  combined 
data  set.  A  computer  program  was  written  to  calculate  the  desired 
probability .  , 
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APPENEI2 
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subroutine  probdif ( nz , z . ny . y . xdb , prcb , iset ) 
c 

c  calculates  probability  that  z  distribution  is  z  dB  less 

c  than  y  distribution 

c 

c  nz  is  number  of  points  from  z  distribution 

c  x( . )  are  points  from  z  distribution 

c  ny  is  number  of  points  from  y  distribution 

c  y( • )  are  points  from  y  distribution 

c  zdb  is  the  db  difference  to  be  tested 

c  prob  is  the  probability  of  a  zdb  difference 

c  iset-1  indicates  that  points  have  already  been  ordered 

c 

dimension  z( 100) .y( 100) ,zo( 100) . yo( 100) 
if  (iset.eq.l)  go  to  70 
n-nz+ny 
c 

c  ordering  of  z  points 

0 

do  5  i-l.az 
3  zo(l)-z(i) 

if  (nz.eq.l)  go  to  30 

do  20  i-l,nz-l 

do  10  J-i+-l.nz 

if  (zo^)  .It.zo(J))  go  to  10 

t-zo( J ) 

zo( j )-zo(i) 

xo(i)-t 

10  continue 

20  continue 

c 

ordering  of  the  y  points 

30  do  35  i-l.ny 

35  Yo(i)-y(i'; 

if  (ny.eq.l)  go  to  70 
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do  60  i=l,nY“l 
do  40  j=i-^I.ny 
if  (yoCi) . It . yoC j ) )  go  to  40 
t=yo( j ) 
yo(J)=yo(i) 
yo(i)=t 
40  Gontinue 
60  continue 

70  continue 
c 

c  nsmallz  is  the  *  of  z  points  less  than  yo(l)-xdb 

c  nsmally  is  the  ♦  of  y  points  less  than  zoCD-rzdb 

c 

nsmallz-O 

nsmally-O 

pl-l.O 

do  80  i-l.nz 

if  (zo(i).gt.yo(l)-xclb)  go  to  85 
nsinallz=-nsinallz+ 1 
80  continue 

85  if  (nsmallz. eq. 0)  go  to  95 

go  to  120 

95  do  100  j-l.ny 

if  (zo(l)  .It.yoC  j)-zcib)  go  to  105 
ns!nally-nsinally+ 1 
100  Continue 

105  continue 

0 

o  nlargey  is  the  *  of  y  points  greater  than  zo(az)-zdb 

c  nlargez  is  the  *  of  z  points  greater  than  yo(ny(-zcib 

c 

120  nlargey-0 

nlargex-0 

do  130  j-l,ny 
Jj-ny+l-j 

if  (yo( j j ) . It . zo(nz)-xdb)  go  to  135 
nlargey-nlargey-^  1 
130  continue 

135  if  C nlargey . eq . 0 )  go  to  145 
go  to  170 
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o  o  o 


145  ,  do  150  i=l . D2 

ii=QZ+l-i 

if  ( yoCny)  • 
nlargex=Qlargex+ 1 
150  continue 

155  continue 


170  continue 


calculation  of  probabillitY 


50 


G 

O 

c 

200 


prob-0 . 0 
nxp-nx+1 
iiYp*nY+l 

if  (nsmallx.gt.O.and.nlargey.gt.O)  go  to  200 
if  (nsmallY-gt .Q.and.nlargex.gt.O)  go  to  300 
if  (nsmallx.gt .O.and.nlargex.gt.O)  go  to  400 
if  ( nsmall Y • gt . o . and . nlargey • gt ■ 0 )  go  to  500 
print  50 

format  ('  «**•«»  impossible  condition  ) 
return 

case  #1  nsmallx>0  and  nlargeY>0 
f-1. 


do  250  i-nsmallx , nx 
if  (i.gt .nsmallx)  go  to  210 

do  205  j-1, nsmallx 
nxp-nxp-1 
np-np-1 
f-(f •nxp)/np 
205  continue 

if  (i.lt.nx)  go  to  210 
208  b-1. 

go  to  246 

210  if  (i.eq.nx)  go  to  208 
nxp-nxp-1 
np-np-1 
f-(f *nxp)/ np 
m-np-1 
nYp-nY 

g-(l .  *nYp)''ffl 

do  240  J-nlargeY.ny-l 

if  ( j -gt .nlargeY)  go  to  230 


II 


o  o  o 


do  225  k=l.nlargeY 

nyp-nyp-l 

m-m-l 

g-(g*nyp)/m 

225  continue 

go  to  240 
230  nyp-nyp-1 
m-m-1 

g-Cg*nYp)/m 

li-h+g 

240  continue 

245  prob-prob+f *h 
250  continue 

prob-1 . -prob 
return 

case  *2  nsmaLlly^O  and  nlargez>0 

300  continue 
f-1. 

do  350  i-nsfflallY.nY 
if  (i.gt .nsoallY)  go  to  310 

do  305  J-l,nsmallY 


• 

nYP-nYP~l 

np-np-1 

f-(f»nYP)/nip 

305 

continue 

if  Ci.lt.nY) 

go  to  310 

308 

b-l. 

go  to  345 

310 

if  (i.eq.nY) 
nYP-nYP-1 
np-np-1 
f-(f *nYP)/np 

go  to  308 

320 

m-np-l 

nzp-nx 

g-(  1 .  *nxp)/in 

do  340  j-nlargex.nx-l 
if  (J .gt .nlargex)  go  to  330 
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do  325  l£=l,nlargex 

nxp=nxp-l 

m-m-1 

g-(g*n2p)/m 

h-g 

325  continue 

go  to  340 
330  nxp-nzp-l 

g-(g*naT))  /m 
h-h+g 

340  continue 

345  prob-prolj+f  *h 

350  continue 

return 

case  #3  nsaialix>0  and  nlargey'O 

400  continue 

if  (nsmallx.ne.nlargex)  go  to  401 

prob-0 . 5 

rettxrn 

401  ns-iain(nsinallx.nlargex) 
nl-niax(  nsmallx .  nlar-  ex ) 
f-1.0 

do  450  i-nl,nx 

if  (1. gt.nl)  go  to  410 

do  405  J-l,nl 
nxp«nxp-l 
np-np- 1 
f-(f *nxp) /np 
405  continue 

if  (i.lt.nx)  go  to  420 
408  d-l. 

go  to  445 
410  nxp«nxp-l 
np-np-1 
f-(f *nxp) / np 
if  (i.eq.ax)  go  to  408 
420  d-O . 

m-np-l 

nyp-ny 

g-(l.*nyp)/m 

if  (nyp.eq.l)  go  to  422 

m-m-l 

nyp-nyp- 1 

g-(g’'nyp)/iD 
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9  m'  9  '  f"  • 


422  f 1=1.0 

nxpp=n2:p 

do  425  i=*l,ns 
nxpp-nxpp'l 

if  (nxpp.eq.O)  go  to  426 
a-m-l 

f l-(f l*nxpp) /m 
h-h+g’fl 

425  continue 

426  m-np-l 
g-(nxp-l.  )/in 
nyp-QY 

m-m-l 

g-(g»nYp)/m 

fl-1. 

do  440  j-l.ny-l 

nyp-nyP"! 

m-m-l 

fl-(fl*nYp)/ffl 
h-ii+g'fl 
440  continue 

445  prob-prob+f  *h. 

450  continue 

if  (nsmallx.gt .nlargex)  prob-1 . O-prob 
return 

0 

c  case  *4  nsmally>0  and  nlargey’O 

c 

500  continue 

if  (nsmally.ne.nlargeY)  go  to  501 

prob-0 . 5 

return 

501  ns-minCnsmallY.nlargeY) 
nl-maxC  nsmallY . nlargey) 
f-1.0 

do  550  i-nl , ny 

if  (i. gt.nl)  go  to  510 

do  505  j -l.nl 
nyp-nyp-l 
np-np-1 
f-(f *nYp)/ap 
505  continue 
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do  540  J-l,nz-l 
nacp-nxp-1 

n-B-i 

fl-(fl*iixp)/n 
h-h+g’fl 
540  continue 

545  prob-prob+f  *li 
550  continue 

if  (nsmallY.lt.nlargey)  prob=l . 0-prob 

return 

end 
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